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Abstract: Single trial analyses of ensemble activity in alert animals demonstrate that cortical circuits dy¬ 
namics evolve through temporal sequences of metastable states. Metastability has been studied for its potential 
role in sensory coding, memory and decision-making. Yet, very little is known about the network mechanisms 
responsible for its genesis. It is often assumed that the onset of state sequences is triggered by an external 
stimulus. Here we show that state sequences can be observed also in the absence of overt sensory stimulation. 
Analysis of multielectrode recordings from the gustatory cortex of alert rats revealed ongoing sequences of 
states, where single neurons spontaneously attain several firing rates across different states. This single neuron 
multi-stability represents a challenge to existing spiking network models, where typically each neuron is at 
most bi-stable. We present a recurrent spiking network model that accounts for both the spontaneous genera¬ 
tion of state sequences and the multi-stability in single neuron firing rates. Each state results from the activation 
of neural clusters with potentiated intra-cluster connections, with the firing rate in each cluster depending on 
the number of active clusters. Simulations show that the models ensemble activity hops among the different 
states, reproducing the ongoing dynamics observed in the data. When probed with external stimuli, the model 
predicts the quenching of single neuron multi-stability into bi-stability and the reduction of trial-by-trial vari¬ 
ability. Both predictions were confirmed in the data. Altogether, these results provide a theoretical framework 
that captures both ongoing and evoked network dynamics in a single mechanistic model. 
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1 Introduction 

Sensory networks respond to stimuli with dynamic and time-varying modulations in neural activity. Single 
trial analyses of sensory responses have shown that cortical networks may quickly transition from one state 
of coordinated firing to another as a stimulus is being processed or a movement is being performed [1-7]. 
Applying a Hidden Markov Model (HMM) analysis to spike trains recorded from neuronal ensembles in the 
gustatory cortex (GC) revealed that taste-evoked activity progresses through a sequence of metastable states 
[8]. Each metastable state can be described as a pattern of ensemble activity (a collection of firing rates) 
lasting tens to hundreds of milliseconds. Transitions between different states are characterized by coordinated 
changes in firing rates across multiple neurons in an ensemble. This description of taste-evoked activity in 
GC has been successful in capturing essential features of gustatory codes, such as trial-to-trial variability [8] 
and learning-dependent changes [9], In addition, describing neural responses as sequences of metastable states 
has been proposed as an explanation for the role of GC in digestive decisions [10]. The observation that 
similar dynamics exist also in frontal [1, 2, 11], motor [12], premotor, and somatosensory [7] cortices further 
emphasizes the generality of this framework. While the functional significance of these transitions among 
metastable states is being actively investigated [9], very little is known about their genesis. A recent attractor- 
based model has successfully demonstrated that external inputs can give rise to transitions of metastable states 
by perturbing the network away from a stable spontaneous state [10]. However, sequences of metastable states 
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could also occur spontaneously, hence reflecting intrinsic network dynamics [13-17]. Intrinsic activity patterns 
have indeed been observed in vivo during periods of ongoing activity, i.e., neural activity in the absence of 
overt sensory stimulation. Such ongoing patterns have been reported in visual [18-21], somatosensory [22, 23] 
and auditory [24] cortices, and the hippocampus [24-26], The functional significance of ongoing activity, its 
origin and its relationship to evoked activity remain, however, elusive. Here, we analyzed ongoing activity to 
investigate whether transitions among metastable states can occur in the absence of sensory stimulation in the 
GC of alert rats. We found that ongoing activity, much like evoked activity, is characterized by sequences of 
metastable states. Transitions among consecutive states are triggered by the co-activation of several neurons in 
the ensemble, with some neurons capable of producing multiple i.e., 3 or more firing rates across the different 
states (a feature hereby referred to as multi-stability). To elucidate how a network could intrinsically generate 
transitions among multi-stable states, we introduce a spiking neuron model of GC capable of multi-stable states 
and exhibiting the same ongoing transition dynamics observed in the data. We prove the existence of multi¬ 
stable states analytically via mean field theory, and the existence of transient dynamics via computer simulations 
of the same model. The model reproduced many properties of the data, including the pattern of correlations 
between single neurons and ensemble activity, the multi-stability of single neurons, its reduction upon stimulus 
presentation, and the stimulus-induced reduction of trial-by-trial variability. To our knowledge, this is the first 
unified and mechanistic model of ongoing and evoked cortical activity. 

2 Materials and Methods 

2.1 Experimental subjects and surgical procedures 

Adult female Long Evans rats (weight: 275 — 350 gr) were used for this study [27]. Briefly, animals were kept in 
a 12h:12h light/dark cycle and received ad lib. access to food and water, unless otherwise mentioned. Movable 
bundles of sixteen microwires (25 fim nichrome wires coated with fromvar) attached to a mini-microdrive 
[27, 28] were implanted in GC (AP 1.4, ML ±5 from bregma, DV 4.5 from dura) and secured in place with 
dental acrylic. After electrode implantation, intra-oral cannulae (IOC) were inserted bilaterally and cemented 
[29, 30], At the end of the surgery a positioning bolt for restraint was cemented in the acrylic cap. Rats were 
given at least 7 days for recovery before starting the behavioral procedures outlined below. All experimental 
procedures were approved by the Institutional Animal Care and Use Committee of Stony Brook University and 
complied with University, state, and federal regulations on the care and use of laboratory animals. More details 
can be found in [27]. 

2.2 Behavioral training 

Upon completion of post-surgical recovery, rats began a water restriction regime in which water was made 
available for 45 minutes a day. Body-weight was maintained at > 85% of pre-surgical weight. Following few 
days of water restriction, rats began habituation to head-restraint and training to self-administer fluids through 
IOCs by pressing a lever. Upon learning to press a lever for water, rats were trained that water-reward was 
delivered exclusively by pressing the lever following an auditory cue (i.e., a 75 dB pure tone at a frequency of 5 
kHz). The interval at which lever-pressing delivered water was progressively increased to 40±3 s (ITI). In order 
to receive fluids rats had to press the lever within a 3 s long window following the cue. Lever-pressing triggered 
the delivery of fluids and stopped the auditory cue. During training and experimental sessions additional tastants 
were automatically delivered at random times near the middle of the ITI, at random trials and in the absence of 
the anticipatory cue. Upon termination of each recording session the electrodes were lowered by at least 150/nn 
so that a new ensemble could be recorded. Stimuli were delivered through a manifold of 4 fine polyimide 
tubes slid into the IOC. A computer-controlled, pressurized, solenoid-based system delivered ~ 40 /il of fluids 
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(opening time ~ 40 ms) directly into the mouth. The following tastants were delivered: 100 mM NaCl, 100 
mM sucrose, 100 mM citric acid, and 1 mM quinine HC1. Water (~ 50 /il) was delivered to rinse the mouth 
clean through a second IOC five seconds after the delivery of each tastant. Each tastant was delivered for at 
least 6 trials in each condition. 


2.3 Electrophysiological recordings 

Single neuron action potentials were amplified, bandpass filtered (at 300 — 8000 Hz), digitized and recorded 
to a computer (Plexon, Dallas, TX). Single units of at least 3:1 signal-to-noise ratio were isolated using a 
template algorithm, cluster cutting techniques and examination of inter-spike interval plots (Offline Sorter, 
Plexon, Dallas, TX). 


2.4 Data analysis 

All data analysis and model simulations were performed using custom software written in Matlab (Mathworks, 
Natick, MA, USA), Mathematica (Wolfram Research, Champaign, IL), and C. Starting from a pool of 299 sin¬ 
gle neurons in 37 sessions, neurons with peak firing rate lower than 1 Hz (defined as silent) were excluded from 
further analysis, as well as neurons with a large peak around the 6 — 10 Hz in the spike power spectrum, which 
were considered somatosensory [27, 31, 32], Only ensembles with more than two simultaneously recorded neu¬ 
rons were included in the rest of the analyses (27 ensembles with 167 neurons). We analyzed ongoing activity 
in the 5 seconds interval preceding either the auditory cue or taste delivery, and evoked activity in the 5 seconds 
interval following taste delivery exclusively in trials without anticipatory cue. 


2.4.1 Hidden Markov Model (HMM) analysis 


The details of the HMM have been reported in e.g. [1, 7, 8, 33]; here we briefly outline the methods of this 
procedure, especially where our methods differ from previous accounts. Under the HMM, a system of N 
recorded neurons is assumed to be in one of a predetermined number of hidden (or latent) states. Each state is 
defined as a vector of N firing rates, one for each simultaneously recorded neuron. In each state, the neurons 
were assumed to discharge as stationary Poisson processes (“Poisson-HMM”). 

We denote by [j/j( 1),..., yi(T)\ the spike train of the z-th neuron in the current trial, where yi(t) = k if k 
spikes occur in the interval [t, t + dt] and zero otherwise; dt = 1 ms and T is the total duration of the trial in 
ms (or, alternatively, the number of bins since dt = 1 ms). Denoting with St the hidden state of the ensemble 
at time t, the probability of having k spikes/s from neuron i in a given state m in the interval dt is given by the 
Poisson distribution: 


p(yi(t) = k\S t = to) 


( v i (m)dt) k e- Vi{ - m ) dt 
k\ 


( 2 - 1 ) 


where i/j(m) is the firing rate of neuron i in state to. The firing rates z/, (to) completely define the states and 
are also called emission probabilities in HMM parlance. We matched the HMM model to the data segmented 
in 1 ms bins. Given the short duration of the bins, we assumed that in each bin either zero or one spikes are 
emitted by each neuron, with probabilities p(yi(t) = 0|St = m) = e~ Vi< ' rn ' , ' dt and p(yi(t) = 1|5* = m) = 
1 — e ~ v i { m )- dt , respectively. We also neglected the simultaneous firing of two or more neurons (a rare event): 
if more than one neuron fired in a given bin, a single spike was randomly assigned to one of the firing neurons. 

The HHM model is fully determined by the emission probabilities defining the states and the transition 
probabilities among those states (assumed to obey the Markov property). The emission and transition prob¬ 
abilities were found by maximization of the log-likelihood of the data given the model via the expectation- 
minimization (EM), or Baum-Welch, algorithm [34]. This is known as “training the HMM.” For each session 
and type of activity (ongoing vs. evoked), ensemble spiking activity from all trials was binned at 1 ms intervals 
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prior to training assuming a fixed number of latent states M [8, 33], We used always the same state as the initial 
state (e.g., state “1”), but ran the algorithm starting 200 ms prior to the interval of interest to obtain an unbiased 
estimate of the first state of each sequence. For each given number of states M, the Baum-Welch algorithm was 
run 5 times, each time with random initial conditions for the transition and emission probabilities. The proce¬ 
dure was repeated for M ranging from 10 to 20 for ongoing activity trials, and from 10 to 40 for evoked activity 
trials, based on previous accounts [8, 10, 33]. For evoked activity, each HMM was trained on all four tastes 
simultaneously, using overall the same number of trials as for ongoing activity. Of the models thus obtained, 
the one with largest total likelihood was taken as the best HMM match to the data, and then used to estimate 
the probability of the states given the model and the observations in each bin of each trial (a procedure known 
as “decoding”). During decoding, only those states with probability exceeding 80% in at least 50 consecutive 
1-ms bins were retained (henceforth denoted simply as “states”). 

For the firing rate analysis of Fig. 2, the firing rates Vi (m) were obtained from the analytical solution of 
the maximization step of the Baum-Welch algorithm. 


Vi(m) 



{ EL r m (t) ) ' 


( 2 . 2 ) 


Here, r m (t) = p (St = m\y(l ),..., y(T)) is the probability that the state at time t is m, given the observations. 
Note that this procedure allows for sampling the variability of Vi(m) across trials, which can be then used for 
inferential statistics (see below, “State specificity and firing rate distributions”). 

The distribution of state durations across sessions, t s , was fitted with an exponential function f(t s ;a, b) = 
a ■ exp (bt s ) using non-linear least squares minimization and 95% confidence intervals (CIs) were computed 
according to a standard procedure (see e.g. [35], p. 50). 


2.4.2 Firing rate modulations and change-point analysis 

We detected changes in single neurons’ instantaneous firing rate (in single trials) by using the change point 
procedure based on the cumulative distribution of spike times of [36]. Briefly, the procedure selects a sequence 
of putative variations of firing rate, called “change points,” as the points where the cumulative record of the 
spike count produces local maximal changes in slope. A putative change point was accepted as valid if the 
spike count rates before and after the change point were significantly different (binomial test, p < 0.05). See 
[36] for details. This method, which is based on a binomial test for the spike count [36], was chosen because 
it is a simple and general tool for analyzing cumulative records, and had been already empirically validated for 
the analysis of change points in GC [37], 

2.4.3 Single neuron responsiveness to stimuli 

Neurons with significant post-stimulus modulations of firing rate were defined as stimulus-responsive. The 
analysis was based on the CP procedure performed in [—1, 5] s intervals around stimulus delivery occurring at 
time zero (see [37] for further details). The CP procedure was performed on the cumulative record of spike 
counts over all evoked trials for a given taste. If any CP was detected before taste delivery, the CP analysis 
was repeated with a more stringent criterion, until no pre-stimulus change points were detected. If no CP was 
found in response to any of the four tastes, the neuron was deemed not responsive. If a CP was detected in the 
post-delivery interval for at least one taste, a t-test was performed to establish whether the post-CP activity was 
significantly different from baseline activity during the 1 s interval before taste delivery (p < 0.05), and only 
then accepted as a significant CP. Neurons with at least one significant CP were considered stimulus-responsive. 
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2.4.4 Change point-triggered average (CPTA) 


The CPTA of ensemble transitions, formally reminiscent of the spike-triggered average [38, 39], estimates the 
average HMM transition at a lag t from a CP. Given a set of CPs at times / , ..... t p \ for the i-th neuron, its 
CPTA is defined as 


CPTA (i) (f) 



E 

p=i 


[t+tp ’ ,t+tp 5 +At] 


(2.3) 


where (...) denotes average across trials. At = 50 ms, and I\ a M = 1 if a transition occurs in the interval [a, 6 ], 
and zero otherwise. Significance of the CPTA was assessed by comparing the simultaneously recorded data with 
100 trial-shuffled ensembles (t-test, p < 0.05, Bonferroni corrected for multiple bins), in which we randomly 
permuted each spike train across trials, in order to scramble the simultaneity of the ensemble recordings. A 
CPTA peak at zero indicates that, on average, CPs co-occur with ensemble transitions, while e.g. a significant 
CPTA for positive lag indicates that CPs tend to precede transitions (see Fig. 2B). 


2.4.5 State specificity and firing rate distributions 

We defined a neuron as “state-specific” if its firing rate varied across different HMM states. In order to assess 
state specificity for the i-th unit, we collected the set of firing rates across all trials and states (see earlier Section 
“Hidden Markov Model (HMM) analysis”). To determine whether the distribution of firing rates varied across 
states, we performed a non-parametric one-way ANOVA (unbalanced Kruskal-Wallis, p < 0.05). The smallest 
number of significantly different firing rates for a given neuron in a given state was found with a post-hoc 
multiple comparison rank analysis (with Bonferroni correction). Given a p-value pik for the pairwise post-hoc 
comparison between states l and fc, we considered the symmetric M x M matrix A with elements Aik = 1 if 
the rates were different (p^ < 0.05) and .4/ fc = 0 otherwise. For example, consider the case of 4 states and the 
following two outcomes of multiple comparisons: 


= 

'■Oil' 

••11 

A< 2 ) = 

'•OOF 

••01 


... 1 


... 0 


(2.4) 


In the first case, the firing rate in state 1 is different than in states 3 and 4 (first row); it is different in state 2 
vs. 3 and 4 (second row), and in state 3 vs. 4 (third row). The conclusion is that /i / /a, f \ / / 4 (from the 
first row), and /3 7 ^ / 4 (from the 3rd row), while / 2 , although different from f:> and / 4 , may not be different 
from fi. It follows that at least 3 different firing rates are found across states. In the example of matrix ,4 ( ' 2 \ 
the smallest number of different rates is 2 (e.g., fi / 4 ). 

2.4.6 Single unit-ensemble correlations and 0 statistics 

We estimated the nonparametric correlation between ensemble and single units firing rate modulations in single 
trials, by means of their (j> statistics. For each ensemble transition between consecutive states rri -0 n, we 
looked at whether the i-th neuron increased (s*(to — 0 n) = + 1 ) or decreased (s*(m -0 n) = — 1 ) its firing 
rate during the transition, and analogously for the ensemble mean firing rate s ens (m —0 n) = ±1. For each 
transitions to —>■ n, we counted the number of ensemble neurons that increased their firing rate when the whole 
ensemble increased its firing rate 


N 

n ++ (m —0 n) = ^ (s*(to -0 n) = +1 & s eas (m -0 n) = +1) , (2.5) 

i=l 
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and analogously for the other three possibilities n.\ _, n _ n _, thus obtaining a 2 x 2 matrix n a b(m —> n), 

for a,b = +, —. We then compiled the contingency matrix with elements N a (> = ^ n a b(m —► n), by summing 
each n a f,(m —> n) over all transitions and sessions. Based on these quantities, the ^-statistics [40, 41] was 
computed as 


_det(iV ab )_ 

\/N+b)(J2b ^o+)(Ea N a -) 


( 2 . 6 ) 


If changes in single neurons perfectly matched ensemble changes (as in the case of UP and DOWN states), then 
(f> = 1; while (j> = 0 in case of no correlation and 0 < (f> < 1 in case of partial co-activation of neurons and 
ensembles. 


2.4.7 Fano factor 

The raw Fano factor (FF) for the z-th unit is defined as 

var[iVj(f)] 

mt)) 


(2.7) 


where A r , ; (t) is the spike count for the z-th neuron in the [t, t + At] bin, and (...,) and var are the mean and 
variance across trials. In order to control for the difference in firing rates at different times, the raw FF was 
mean-matched across all bins and conditions, and its mean and 95% Cl were estimated via linear regression 
(see [42] for details). We used a window of 200 ms for the comparison between model and data. Results were 
qualitatively similar for a wide range of windows sizes (At = 50 to 500 ms), sliding along at 50 ms steps. 


2.5 Spiking neuron network model 

Our recurrent network model comprised N = 5000 randomly connected leaky integrate-and-fire (LIF) neurons, 
with a fraction he = 0.8 of excitatory (E) and nj = 0.2 inhibitory (I) neurons. Connection probability pp a 
from neurons in population a £ E, I to neurons in population (3 £ E,I were pee = 0-2 and pei = Pie = 
pn = 0.5. A fraction / = 0.9 of excitatory neurons were arranged into Q different clusters, with the remaining 
fraction belonging to an unstructured (“background”) population [43], Synaptic weights Jp a from neurons 
in population a £ E,I to neurons in population (3 £ E, I scaled with N as Jp a = jp a /VN, with jg a 
constants drawn from a normal distribution with the following mean values (units of mV): jEi = 3.18, jiE = 
1.06, jn = 4.24, Jee = 1.77; and variance Var[J] = S 2 J 2 , with S 2 = 0.01. Within an excitatory cluster 
synaptic weights were potentiated, i.e. they took values J+Jee with J + > 1 , while synaptic weights between 
neurons belonging to different clusters were depressed to values J-Jee, with J_ = 1 — 7 (J+ — 1 )f/Q < 1 , 
with 7 = 0.5. The latter relationship between J + and ./_ insures balance between overall potentiation and 
depression in the network [43]. Below spike threshold, the membrane potential V of each LIF neuron evolved 
according to 

An j, ~ C T T m (/ rec T / ex t -f fstim) > (2-8) 

at 

with a membrane time constant r m = 20 ms for excitatory and 10 ms for inhibitory neurons. The input current 
was the sum of a recurrent input / rec , an external current / exl representing an ongoing afferent input from other 
areas, and an input stimulus J st i m representing a delivered taste during evoked activity only. In our units, a 
membrane capacitance of 1 nF is set to 1. A spike was said to be emitted when V crossed a threshold V[hr. after 
which V was reset to a potential V[ e set = 0 for a refractory period of r re f = 5 ms. Spike thresholds were chosen 
so that, in the unstructured network (i.e., with J + = </_ = 1), the E and I populations had average firing rates 
of 3 and 5 spikes/s, respectively [43], The recurrent synaptic input /,% to neuron i evolved according to the 
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dynamical equation 


rip 

u,± rp 


dt 


N 

E 

3 =1 


= -4c + E^E^-4) 

k 


(2.9) 


where tj, was the arrival time of fc-th spike from the j-th pre-synaptic neuron, and t s was the synaptic time 
constant (3 and 2 ms for E and I neurons, respectively). The PSC elicited by a single incoming spike was 
-±i_ e ~ t / T ‘Q{t), where 0(f) = 1 for t > 0, and 0(f) = 0 otherwise. The ongoing external current to a neuron 
in population a was constant and given by 


(cxl — 'A>0^ext - (2.10) 

where N ext = n E N , p a0 = p EE , J a o = jao/VN with j E0 = 0.3, j I0 = 0.1, and z/ ext = 7 spikes/s. During 
evoked activity, stimulus-selective neurons received an additional transient input representing one of the four 
incoming stimuli. The percentage of neurons responsive to one, two, three or four stimuli was modeled after 
the estimates obtained from the data, which implied that stimuli targeted overlapping neurons (see Results). We 
tested two alternative model stimuli: a biologically realistic stimulus z/* im (f) resembling thalamic stimulation 
[44], modeled as a double exponential with peak amplitude of 0.3zz ex t and rise times of 50 ms and decay times 
of 500 ms, or a stimulus of constant amplitude ranging from 0 to 0.5zz ex t (“box” stimulus). In the following 
we measure the stimulus amplitude as percentage of z/ ext (e.g., “30%” corresponds to z/ stim = 0.3z/ cxt . The onset 
of each stimulus was always f = 0, the time of taste delivery. The stimulus current to a neuron in population a 
was constant and given by I stim = N ext p a0 J a0 p stim . 

2.5.1 Mean field analysis of the model 

The spiking network model described in the previous subsection is a complex system capable of many behav¬ 
iors depending on the parameter values. One main aim of the model is finding under what conditions it can 
sustain multiple configurations of activity that can be later interpreted as HMM states. Parameter search was 
used relying on an analytical procedure for networks of LIF neurons known as “mean field theory” or “pop¬ 
ulation density approach” (see e.g. [43, 45, 46]). Under the conditions stated below, this theory provides a 
global picture of network behavior together with the associated parameter values, which can then be tested in 
model simulations. Under typical conditions, each neuron of the network receives a large number of small 
post-synaptic currents (PSCs) per integration time constant. In such a case, the dynamics of the network can be 
analyzed under the diffusion approximation, which is amenable to the population density approach. The net¬ 
work has a = 1,, Q + 2 sub-populations, where the first Q indices label the Q excitatory clusters, a = 0+1 
labels the “background” excitatory population, and a = Q+2 labels the homogeneous inhibitory population. In 
the diffusion approximation [47-49], the input to each neuron is completely characterized by the infinitesimal 
mean //„ and variance cr ( 2 l: of the post-synaptic potential. Adding up the contributions from all afferent inputs, 
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/i„ and cr„ for an excitatory neuron belonging to cluster a are given by [43] 


fJ'Ct — TV 




Q -1 


PEEJ+jEE^a + ^ PEEJ-jEEVp I + TIe{\ — f)pEEJ-jEEV 


(bg) 

E 


— niPEljElVl + TlEPEojEOVext 


)■ 


^Q; - T" m 


-(¥ (’ 


Q-i 


p_E£)(J + j_E_E) 2 (l + S 2 )l/ a + ^ P_E_e(J_Jbb) 2 (1 + (5 2 )^ 


+ri£)(l - f)pEE{J-jEE) 2 { 1 + & 2 )ve 9) + n I p E ij 2 EI {l + <5 2 w) , 


where is the firing rate of the unstructured (“background”) E population. Afferent current and variance to 
a neuron belonging to the background E population and to the homogeneous inhibitory population are 


r Q 

0 E ^ = r m,.EV A /vf—^ PEEJ-jEEVp + n E { 1 — f)PEEjEEV E 3 ^ 

^ 0=1 

~niPEljEI^I + UEPEojEOVext'j , 

{v E 9) Y = T rn>E C % ^f'52pEE{J-jEE) 2 ( 1 + 5 2 )V0 + n E { 1 - f)pEE3 2 E E ( 1 + ^) v E^ 

^ 0=1 

+niPEij 2 Ei( 1 + <^Vj) , 

01 = T m jVN y ^PlEjlEb'0 +71_e(1 - f)PlEjlE^Y : 9 ' > 

^ 0=1 

-nipujuvi + nsp/oj/o^ext) , 

cr? = T mi /( n ^y^p /£ ;j^(l + ( 5 2 )l/ / 3 + 71^(1 -/)p 7 E jf E (l + 5 2 )^ ff) 

^ / 3=1 

+rc/P//j 2 /(l + <5 2 )^/)- (2.11) 

Parameters were chosen so as to have a balanced unstructured network. In other words, our network with 
J + = J_ = 1 (where all E—>E synaptic weights are equal) would operate in the balanced asynchronous 
regime where incoming contributions from excitatory and inhibitory inputs balance out and the peri-stimulus 
time interval over time is approximately flat. In such a regime, one can solve for the excitatory and inhibitory 
firing rates as linear functions of the external firing rate v ex t, up to terms of Of 1 /\/N), provided the connection 
probabilities and synaptic weights satisfy one of the following conditions [50-52]: either 


or 


PElPlEjEljlE 

PEEPlljEE 


PlIjEojll 
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Asynchronous activity was confirmed by computing the network synchrony, defined as [53] 


E = 



(2.14) 


where af is the time-averaged variance in the instantaneous firing rate Vi(t) of the Mil neuron, [...] denotes the 
average over all neurons in the network, and <j^ op is the time-averaged variance of the instantaneous population 
firing rate [vi(t)\. The network synchrony E is expected to be of 0(1) in a fully synchronized network, but of 
order 0(l/y/N) in networks of size N that are in the asynchronous state. We used bins of variable size (from 
10 ms to 50 ms) to estimate E, obtaining similar results in all cases. The unstructured network has only one 
dynamical state, i.e., a stationary point of activity where all E and I neurons have constant firing rate ve and ///, 
respectively. In the structured network (where J + > 1), the network undergoes continuous transitions among 
a repertoire of states, as shown in the main text. Not to confuse the networks states of activity with the HMM 
states, we shall from now on use the term “network configurations” instead. Admissible network configurations 
must satisfy the Q + 2 self-consistent mean field equations [43] 


^a Fa ^aC^)) > 


(2.15) 


where it = \v\, ..., uq, v^ 9 \ Vi] is the firing rate vector and F a (^i a , al) is the current-to-rate response 
function of the LIF neurons. For fast synaptic times, i.e. r s /T m << 1 , F n ( fi lt . a'l ) is well approximated by 
[54, 55] 

/ /-Oefrc* „ \ 1 


F'ot (/^o: ? ^cn) — I 7"ref H - 


' #eff,c 


e“” [1 + erf(u)] 


(2.16) 


where 


(_) _ Vthr ’ a L la l „7. 

\ Llhjcx , 




O OL 

Preset,a /^o 

O nt 


ak n 


where k a = \/r Si0 /r my , is the square root of the ratio of synaptic time constant to membrane time constant, 
and a = |C(l/2)|/-\/2 ~ 1.03. This theoretical response function has been fitted successfully to the firing rate 
of neocortical neurons in the presence of in vivo-like fluctuations [56-58]. The fixed points of the mean 
field equations were found with Newtons method [59], The fixed points can be either stable (attractors) or 
unstable depending on the eigenvalues A Q of the stability matrix 


Fa /3 


1 / dF a ^ a (^),al(lt)) 

T s ,a y dvp 


dF a (^),a 2 a (^)) dal 
dal dvp 



(2.17) 


evaluated at the fixed point 17** [60]. 1 If all eigenvalues have negative real part, the fixed point is stable (attrac¬ 
tor). If at least one eigenvalue has positive real part, the fixed point is unstable. Stability is meant with respect 

^qual indices are not summed over. The published version of the current article in [61] has a typo that was corrected in 2.17. 
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to an approximate linearized dynamics of the mean and variance of the input current: 


dm a 


Ts,a ds 2 a 

2 dt 


v a {t) 


—m a + n a ( 1?) , 

“4 + 4C^) . 


where /;„ and cr„ are the stationary values for fixed nu given earlier. For fast synaptic dynamics in the asyn¬ 
chronous balanced regime, these rate dynamics are in very good agreement with simulations ([62] see [63, 64] 
for more detailed discussions). 


2.5.2 Model simulations 

The dynamical equations of the LIF neurons were integrated with the Euler algorithm with a time step of 
dt = 0.1 ms. We simulated thirty ensembles during ongoing activity and twenty ensembles during evoked 
activity. We chose four different stimuli per session during evoked activity (as in the data). Trials were 5 
seconds long. In all cases, we pre-processed each ensemble by picking at random a distribution of neurons per 
ensemble that matched the distribution of ensemble sizes found in the data, and performed the same analyses 
performed on the data (i.e., HMM, CPTA, FF, and so on). The range of hidden states M for the HMM in the 
model was the same as the one used for the data (see above), except in the case of analyses involving simulated 
ensembles of 30 neurons (Fig. 5A), where M varied from 20 to 60 states. To determine the number of active 
clusters at any given time (Fig. 4), a cluster was considered active if its population firing rate was higher than 20 
spikes/s in a 50 ms bin. This criterion gave a good separation of the firing rates of active and inactive clusters in 
the multi-stable region (Fig. 3B, region with J + > ./'). A cluster was active in a particular bin if its population 
firing rate crossed the 20 spks/s threshold in that bin. 


3 Results 

3.1 Characterization of ongoing cortical activity as a sequence of ensemble states 

We analyzed ensemble activity from extracellular recordings in the gustatory cortex (GC) of behaving rats. In 
our experimental paradigm, rats were delivered one out of four tastants (sucrose, sodium chloride, citric acid and 
quinine, denoted respectively as S, N, C, and Q) through an intra-oral cannula, followed by a water rinse to clean 
the mouth from taste residues [27]. In half of the trials, taste delivery was preceded by an auditory cue (75 db, 
5 kHz), signaling taste availability upon lever press; in the other half of the trials, taste was delivered randomly 
in the absence of an anticipatory cue. Expected and unexpected trials were randomly interleaved. Rats were 
engaged in a task to prevent periods of rest. We analyzed ongoing activity in the 5 seconds interval preceding 
either the auditory cue or taste delivery, and evoked activity in the 5 seconds interval following taste delivery 
in trials without anticipatory cue. The interval length was chosen based on the fact that population activity in 
GC encodes taste-related information for at least five seconds after taste delivery [37] and the ongoing interval 
was designed to match the length of the evoked one. Visual inspection of ensemble spiking activity during 
ongoing inter-trial periods (Fig. 1) reveals that several neurons simultaneously change their firing rates. This 
suggests that, similarly to taste-evoked activity [8], ongoing activity in GC could be characterized in terms of 
sequences of states, where each state is defined as a collection of constant firing rates across simultaneously 
recorded neurons (Fig. 1; see Appendix 2 for details). We assumed that the dynamics of state transitions 
occurred according to a Markov stochastic process and performed a hidden Markov model (HMM) analysis in 
single trials [2, 7, 8, 33]. For each neural ensemble and type of activity (ongoing vs. evoked), spike trains were 
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Figure 1. Ongoing activity is characterized by sequences of states. Four representative trials from one ensemble of nine 
simultaneously recorded neurons, segmented according to their ensemble states (HMM analysis). Thin black vertical lines 
are action potentials. States are color-coded; smooth colored lines represent the probability for each state; shaded colored 
areas indicate intervals where the probability of a state exceeds 80%. Below each single-trial population raster, average 
firing rates across simultaneously recorded neurons are plotted (states are color-coded). X-axis for population rasters: time 
preceding the next event at (0 = stimulus delivery); X-axis for average firing rates panels: firing rates (spks/s); Y-axis for 
population rasters: left, ensemble neuron index, right, probability of HMM states; Y-axis for firing rate panels: ensemble 
neuron index. 


fitted to several Poisson-HMMs that differed for number of latent states and the initial conditions. The model 
providing the largest likelihood was selected as the best model, but only states with probability larger than 80% 
in at least 50 consecutive 1 ms bins were retained as valid states (called simply “states” in the following; see 
Appendix 2 for details). The number of states across sessions ranged from 3 to 7 (mean±SD: 4.8 ±1.1, Fig. 
2A) with approximately exponentially distributed state durations with mean of 717 ms (95% Cl: [687, 747], 
Fig. 2B). The number of states was correlated with the number of neurons in each state {B 2 = 0.3, p < 0.01). 
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Figure 2. Features of ongoing activity states. A: distribution of number of states per ensemble across all sessions (only 
states lasting longer than 50 ms are included). X-axis: number of states per ensemble; Y -axis: frequency of occurrence. B: 
distribution of state durations across all sessions together with the least squares exponential fit (red line; see Appendix 2). 
X-axis: state duration; X-axis: frequency of states occurrence. 


However, such correlation disappeared when rarely occurring states (occurring for less than 3% of the total 
session duration) were excluded ( R 2 = 0.1, p > 0.1), suggesting that the most frequent states are present even 
in small ensembles. 

To gain further insight into the structure revealed by HMM analysis, we pitted the ensemble state tran¬ 
sitions against modulations of firing rate in single neurons. The latter were found using a change point (CP) 
procedure applied to the cumulative distribution of spike counts [36, 37], Examples from three representative 
trials are shown in Fig. 3A, where it can be seen that CPs from multiple neurons (red dots) tend to align with 
ensemble state transitions (dashed vertical green lines). The relationship between CPs and state transitions was 
characterized with a CP-triggered average (CPTA) of state transitions. The CPTA estimates the likelihood of 
observing a state transition at time t, given a CP in any neuron of the ensemble at time t = 0. If detection of 
CPs and state transitions were instantaneous, a significant peak of the CPTA at positive times would indicate 
that state transitions tend to follow a CP, whereas a peak at negative times would indicate that state transitions 
tend to precede CPs. We found a significant positive CPTA in the interval between t = 0 and t = 200 ms (Fig. 
3B, black trace), with a peak at time t = 0. However, we must be cautious to conclude that state transitions 
tend to follow CPs in this case. Because of our requirement that the states posterior probability must exceed 
80% to detect an ensemble transition, the scored time of occurrence is likely to lag behind the correct time, 
which would skew the CPTA towards positive values. In any case, the CPTA peak is found around t = 0 (Fig. 
3B) and shows that the largest proportion of state transitions co-occurs with a CP. Significance was established 
by comparison with a CPTA in a trial-shuffled dataset (p < 0.05, t-test with Bonferroni correction; see Ap¬ 
pendix 2). In the trial-shuffled dataset, the CPTA was flat indicating that the relation between CPs and ensemble 
transitions was lost (Fig. 3B, blue trace). The fraction of neurons having CPs co-occurring with a state transi¬ 
tion ranged from a single neuron to half of the ensemble (Fig. 3C), indicating that more than one neuron, but 
not all neurons in the ensemble, co-activate during a state transition. In support of this result, we found that 
population firing rate modulations and single neurons modulations were partially correlated across transitions 
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Figure 3. Single neurons co-activation drives ensemble transitions. A: change points (CPs, red dots) and state transitions 
(vertical dashed green lines) are superimposed to a raster plot of simultaneously recorded neurons (same trials as the first 
three in Fig. 1). A'-axis for upper panel: time preceding taste delivery (0 = stimulus delivery); Y-axis: ensemble neuron 
index. B: CP triggered average (CPTA) of ensemble transitions for the empirical dataset (black) and for the trial-shuffled 
dataset (blue). Shaded bounds: standard errors. Thick black line below traces: p < 0.05 (bin-wise t-test with multiple-bin 
Bonferroni correction). A-axis: time lag from CP (s); Y-axis: probability of a transition given a CP at t = 0. C: fraction 
of transitions co-occuring with CPs within a 200 ms window prior to the state transition (significant window from panel B). 
A-axis: number of single neurons CPs divided by the ensemble size; Y-axis: fraction of transitions co-occurring with CPs 
(%). D: correlation (^-statistics) between the signs of single neurons and ensemble firing rate changes for the empirical 
(“Data”) and simulated datasets. In the latter, 90%, 50% and 0% of simulated neurons had firing rate changes matched 
to those of the whole ensemble (boxes in box-plots represent 95% CIs). Y-axis: (^-statistics. E: firing rate distributions 
across states for representative neurons 2 and 3 from the ensemble in Fig. 1 and panel A above (color-coded as in Fig. 1, 
each curve represents the empirical distribution of firing rates in each state). A'-axis: firing rate (spks/s); Y-axis: probability 
density of states. F: number of different firing rates per neuron across all sessions and neurons: 42% of all neurons had 3 
or more different firing rates across states. A-axis: minimal number of different firing rates across states; Y-axis: fraction 
of neurons (%). 
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(Fig. 3D, cj) 2 = 0.21 ± 0.04, mean±SD). We explained the observed correlation by comparing it to three 
surrogate datasets (correlations in the four datasets were all significantly different; X 2 (3) = 110, p < 10- 10 , 
Rruskal-Wallis one-way ANOVA with post-hoc multiple comparisons). The empirical correlation was signif¬ 
icantly larger than in the case of a surrogate dataset in which all neurons changed their bring rates at random 
(cj) 2 = 0.01 ± 0.01, but signibcantly lower than in surrogate datasets with high levels of global coordination. A 
dataset with 90% and 50% of the neurons having co-occurring changes in bring rates yielded a </> 2 = 0.93±0.01 
and <j) 2 = 0.66 ± 0.02, respectively (Fig. 3D). Taken together, these results conhrm that state transitions are 
due to a partial (and variable in size) co-activation of a fraction of neurons in the whole ensemble. Beyond 
unveiling the multi-state structure described above, inspection of the representative examples in Fig. 1 also pro¬ 
vides indications on the spiking behavior of single neurons. Single neurons bring rates appear to be bi-stable in 
some cases and multi-stable in others. To quantify this phenomenon, we computed the bring rate distributions 
across states for each neuron. Fig. 3E shows two representative examples featuring mixtures of distributions 
peaked at multiple characteristic bring rates (different states are color-coded). To bnd out the minimal number 
of bring rates across states for all recorded neurons, we conducted a conservative post-hoc comparison after a 
non-parametric one-way ANOVA of the bring rates for each neuron (see Appendix 2). The mean bring rates 
varied signibcantly across states for 90% of neurons (150/167; Kruskal-Wallis one-way ANOVA, p < 0.05), 
with at least 42% (70/167) of all neurons being multi-stable, i.e., having three or more signibcantly different 
bring rate distributions across states (Fig. 3F). Altogether, these results demonstrate that ongoing activity in 
GC can be described as a sequence of multiple states in which a portion of neurons changes bring activity in 
a coordinated manner. In addition, our analyses show that almost half of the neurons switch between at least 
three different bring rates across states. 

3.2 A spiking network model with a landscape of multi-stable attractors 

We have shown that ongoing activity in rat gustatory cortex is structured and characterized by sequences of 
metastable states. Biologically plausible models capable of generating such states, based on populations of 
spiking neurons, have been put forward [13, 15, 17]. However, the hallmark of these models is bi-stability in 
single neurons, which bre at either low or high bring rate, depending on the state of the network [43, 65, 66 ]. 
Instead, a substantial fraction of GC neurons can exhibit 3 or more stable bring rates (Fig. 3F), raising a 
challenge to existing models. To address this issue, we developed a spiking network model capable of multi¬ 
stability, where each single neuron can attain several bring rates (i.e. 3 or more). The network contains 4000 
excitatory (E) and 1000 inhibitory (I) leaky integrate-and-hre (LIF) neurons (Fig. 4A and Appendix 2) that are 
mutually and randomly connected. A fraction of E neurons form an unstructured (“background”) population 
with mean synaptic strength Jee- The majority of E neurons are organized into Q recurrent clusters. Synaptic 
weights within each cluster are potentiated (with mean value of (J) = J + Jee- with J + > 1), whereas synapses 
between E neurons belonging to different clusters are depressed (mean value (J) = J_Jee , with J_ < 1 ). 
Inhibitory neurons are randomly interconnected to themselves and to the E neurons. All neurons also receive 
constant external input representing contributions from other brain regions and/or an applied stimulus (such as 
taste delivery to mimic our data). As explained in detail below, a theoretical analysis of this model using mean 
held theory predicts the presence of a vast landscape of attractors, and the existence of a multi-stable regime 
where each neuron can bre at several different bring rates. 

3.3 Mean field analysis of the model 

A mean field theory analysis (see Appendix 2) predicts that, depending on the intra-cluster potentiation param¬ 
eter J + , the cortical network can exhibit configurations with one or more active clusters, where we use the term 
“configuration” of the network to distinguish its admissible states of firing rate activities across neurons from 
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Figure 4. Attractor landscape of the spiking network model. A: network architecture. Excitatory (E, triangles) and in¬ 
hibitory (I, squares) LIF neurons are recurrently connected. E neurons organized in Q clusters; intra-cluster synaptic 
connections are potentiated (darker disks), inter-cluster connections are depressed (lighter disks). B: mean held theory anal¬ 
ysis of the network in panel A. Firing rates for excitatory clusters attractor states (diamonds) are shown as a function of the 
intra-cluster potentiation parameter J+ (in units of Jee). Below the critical point J' = 4.2 (i.e. left most vertical dotted 
line) the only stable state has low bring rate ~ 5 spikes/s (blue diamonds). At J+ = J', a bifurcation occurs whereby 
states with active clusters at higher firing rate coexist (upper diamond). There are in total Q such states. As J+ is further 
increased, states with 2 or more active clusters exist. For each state, all active clusters have the same bring rate (reported 
on the vertical axis). The vertical green line at ,7+ = 5.2 represents the value chosen for the spiking network simulations; 
vertical red lines represent critical points, where a new conhguration appears with an increased number of active clusters. 
X-axis: intra-cluster synaptic potentiation (J+); Y-axis: bring rate (spks/s). C: number of active clusters as a function of 
J+ (notations as in panel B). For J+ = J' all attractor states have only one active cluster; for J+ = 5.2 (green line) there 
are 7 possible conbgurations of attractor states, with 1, 2,..., 7 active clusters, respectively. A'-axis: intra-cluster synaptic 
potentiation (J+); Y-axis: bring rate (spks/s). 


the ensemble states revealed by an HMM analysis. A network configuration is fully specified by the number 
and identity of its active clusters, although we shall often focus on the number of active clusters regardless of 
their identity. The predicted population firing rates, as a function of J + and the number of active clusters in 
stable network configurations, are shown in Fig. 4B. For very weak values of J + > 1, the network has one 
stable configuration of activity (“background” configuration), where all E neurons fire at a low rate vl ~ 5 Hz 
(Fig. 4B, blue diamonds left to the leftmost vertical red line at J' = 4.2). As J + Increases beyond the first 
critical point J' > 1, a bifurcation occurs (Fig. 4B). At the bifurcation, a new set of attractors emerges char¬ 
acterized by a single cluster active at a higher firing rate vh , while the rest of the network fires at « zzj, (gray 
diamonds between the two leftmost vertical red lines in Fig. 4B). There are Q possible such configurations, one 
for each cluster. As we increase J + in small steps, we cross several bifurcation points (vertical red lines in Fig. 
4B). To the right of each line, a new set of global configurations is added to the previous ones, characterized 
by an increasing number of simultaneously active clusters, whose firing rates depend on the number of active 
clusters. For example, for J + = 5.2 (green line), the network can be in one of 30 configurations with only 
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one active cluster at v^ 1 = 64 spikes/s; with 2 active clusters at v^ 1 = 62 spikes/s (with ( 3 2 °) = 435 possible 
configurations); with 3 active clusters at = 58 spikes/s (with ( 3 3 °) = 4060 possible configurations), and so 
on. A configuration with q out of Q simultaneously active clusters can be realized in any one of (^) = q qq'^ q y 

possible configurations, giving rise to a vast landscape of attractors. The firing rates u^\q < Q) depend on 
the number q of simultaneously active clusters in decreasing fashion, i.e. v ))' > > ... > > v L . The 

lower firing rate per cluster in a larger number of active clusters is the consequence of recurrent inhibition, since 
the whole network is kept in a state of global balance of excitation and inhibition [50, 51, 65]. For each maxi¬ 
mal number t/ max of allowed active clusters in the range of J + in Fig. 4B, any subset of g max clusters could in 
principle be active, and the network can be in any one of the global configurations characterized by a number of 
active clusters compatible with the intra-cluster potentiation J + . The number of different firing rates per neuron 
in the stable configurations (diamonds) are plotted as a function of J + in Fig. 4C. It is important to realize that 
these configurations are stable, and would thus persist indefinitely, only in a network with an infinite number of 
neurons. Moreover, configurations with the same number of active clusters are equally likely (for a given value 
of J_|_) only if all clusters contain exactly the same number of neurons. We show next in model simulations 
that, in a network with a finite number of neurons, the persistent configurations are spontaneously destabilized, 
resulting in the network going through a sequence of metastable states, as observed in the data. Moreover, slight 
variations in the number of neurons per cluster greatly decrease the number of observed configurations, also in 
keeping with the data. 
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Figure 5. Dynamics of the spiking network model during ongoing activity. Simulation of the network in Fig. 4 with 
4000 excitatory and 1000 inhibitory LIF neurons, Q = 30 clusters at intra-cluster synaptic potentiation J+ = 5.2. A: 
incoming PSC to an excitatory (PSCE, upper plot) and an inhibitory (PSCI, lower plot) neuron: EPSC (blue trace), 1PSC 
(red trace), external current (green line), and total current (black trace), are in a balanced regime. X-axis: time (s); Y - 
axis: post-synaptic current (nA). B: membrane potential from representative excitatory ( Ve , upper plot) and inhibitory (Vj, 
lower plot) neurons (vertical bars: spikes, horizontal dashed lines: threshold for spike emission). X-axis: time (s); Y - 
axis: membrane potential V (mV). For illustration purposes, V was linearly transformed to obtain the threshold for spike 
emission at —45 mV and the reset potential after a spike at —60 mV. 
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3.4 Model simulations: asynchronous activity and multi-stable regime in the spiking network model 

We tested the predictions of the theory (Fig. 4) in a simulation of the model in the regime corresponding 
to J + = 5.2 (green line in Fig. 4B). Ninety percent of the E neurons were assigned to Q = 30 clusters 
with potentiated intra-cluster connectivity J + Jee- The clusters had a variable size with a mean number of 120 
neurons and 1% standard deviation. In the “background” configuration, the network was in the balanced regime 
as shown in Fig. 5A for one E (top) and one I neuron (bottom), respectively. In this regime, the excitatory (blue), 
inhibitory (red), and external (green) post-synaptic currents sum up to a total current with mean close to zero 
(black trace) and large variance. This E/I balance results in irregular spike trains originating from substantial 
membrane potential fluctuations (Fig. 5B). In line with this, the network had a global synchrony index of 
E « 0.01 ~ 1 /vjV [53], which is the signature of asynchronous activity (see Appendix 2). The fluctuations 
generated in the network turned the stable configurations predicted by the mean field analysis (Fig. 4B) into 
configurations of finite duration, and the network dynamically sampled different metastable configurations with 
different patterns of activations across clusters (Fig. 6A), reminiscent of the sequences shown in Fig. 1. In each 
configuration, q out of Q clusters were simultaneously active at firing rate ■ In the example of Fig. 6A, the 
number of active clusters ranged from 3 to 7 at any give time (Fig. 6B). The number of active clusters across all 
sessions was approximately bell-shaped around 4.8 ± 0.9 (meaniSD; Fig. 6C, inset); the configurations with 
appreciable frequency of occurrences had 3 to 7 active clusters, with firing rates ranging from a few to about 
70 spikes/s. Configurations with 1 or 2 active clusters, although predicted in mean field (Fig. 4B), were not 
observed during simulations, presumably due to low probability of occurrence and/or inhomogeneity in cluster 
size (more on this point below). Moreover, the average population firing rate in each cluster was inversely 
proportional to the number of active clusters, as predicted by the theory (Fig. 6C). The actual firing rates as a 
function of the number of active clusters were also in good agreement with the values predicted by mean field 
theory (Fig. 4B). As the network hops across configurations, the firing rates of single neurons also change over 
time, jumping to values that depend on the number of co-active clusters at any given time (Fig. 6D; different 
clusters are in different colors). These results suggest that this model can provide an explanation of the ongoing 
dynamics observed in the data, as we show next. 

3.5 Ensemble states in the network model 

The model configurations with multiple active clusters could be the substrate for the states observed in the 
data. To show this, we performed an HMM analysis on 30 simulated model sessions, each containing 40 
ongoing trials per session, each trial lasting 5 seconds (Fig. 7). Since all neurons in the same cluster tend 
to be active or inactive at the same time (Fig. 6A), it is sufficient to consider one neuron per cluster in the 
HMM analysis. Fig. 7 shows examples of segmentation in states for a simulated ensemble of 30 neurons (each 
randomly chosen from a different cluster). The duration and distribution of states in the sequences depended 
on a number of factors, among which the amount of heterogeneity in cluster size (larger clusters were more 
likely to be in an active state), the synaptic time constants and the intra-cluster potentiation parameter J + . For 
our set of parameters we found 34 ± 4 states per session, ranging from 29 to 39. Examples of these states 
are shown in Fig. 7. Note that different configurations with the same number, but different identity, of active 
clusters would produce different states, thus in principle our network could produce a very large number of 
states if the clusters would transition independently. Three factors limit this number to the actually observed 
one: 1) due to the networks organization in clusters competing via recurrent inhibition, the number of predicted 
configurations that are stable have between 1 and 7 active clusters for J + = 5.2, as predicted by the theory 
(Fig. 4B, diamonds on the green vertical line). No configurations with 8 ,..., 30 clusters are stable, because 
activation of one cluster will tend to suppress, via recurrent inhibition, the activation of other clusters; 2) not 
all predicted configurations have the same likelihood of being visited by the dynamics of the network; in our 
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Figure 6. Dynamics of the spiking network model during ongoing activity. A: a representative rasterplot from excitatory 
clustered neurons (each dot represents a spike, background population not shown). Clusters of neurons that are currently 
active appear as darker regions of the raster. A'-axis: time (s); K-axis: neuron index. B: time course of the number of active 
clusters from the representative trial in panel A. A-axis: time (s); A-axis: number of active clusters. C: average firing 
rates in the active clusters as a function of the number of active clusters across all simulated sessions (bars: S.D.). A'-axis: 
number of active clusters; Y-axis: average cluster firing rate (spks/s). Inset: occurrence of states with different counts of 
active clusters for 5% stimulus amplitude. A'-axis: number of active clusters; Y-axis: frequency of occurrence (% of total 
time). D: instantaneous cluster firing rate in three representative clusters (red, blue, and green lines) from trial in panel A. 
A-axis: time (s); Y-axis: firing rate (spks/s). 


example, the probability of observing configurations with less than 3 active clusters in model simulations have 
very low probability (Fig. 6 C), further reducing the total number of visited configurations; 3) due to a slight 
heterogeneity in cluster size ( 1 % variance across clusters), configurations with the larger clusters active tend to 
be sampled the most (i.e., they tend to recur more often and last longer), presumably due to a larger basin of 
attraction and/or a deeper potential well in the attractor landscape. All these factors working together produce an 
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Figure 7. Ongoing activity in the model is characterized by sequences of states. Two representative trials from an ensemble 
with 30 neurons from Fig. 6A, one neuron per cluster (notations as in Fig. 1 A). Raster plots (upper panel) are overlaid with 
HMM states, together with their firing rates across neurons (lower panel). Upper panel: X-axis, time preceding the next 
event (0 = stimulus delivery); Y-axis, left: ensemble neuron index; right: probability of HMM states. Lower panel: X-axis, 
firing rates (spks/s); Y-axis: ensemble neuron index. 


itinerant network dynamics among only a handful of configurations (compared with the large number possible 
a priori), that in turn originate only a handful of ensemble states in an HMM analysis of the same simulated 
data (Fig. 7). The presence of recurrent inhibition and the heterogeneity in cluster size are also responsible for 
a higher degree of cluster co-activation than expected by chance had the clusters been independent. This was 
tested by comparison with a surrogate dataset obtained by random trial shuffling of the simulated data. The 
distributions of the number of cluster co-activations in 50 ms bins in the original vs. the shuffled data were 
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significantly different (Kolmogorov-Smirnov test, N e g- = 3 x 10 4 , p < 10~ 12 [59]). In particular, activations 
of single clusters were more frequent in the trial-shuffled dataset whereas co-activations of multiple clusters 
were more frequent in the original model simulations (one-way ANOVA with Bonferroni corrected multiple 
comparisons, p < 0.01). Finally, the average number of different configurations was significantly larger in the 
trial-shuffled dataset (Wilcoxon rank-sum test across sessions, z = —2.4, p < 0.05). Taken together, these 
results show that the model network, although producing asynchronous activity hence low pairwise spike train 
correlations (not shown), induces more coordinated clusters co-activations and a much more limited number of 
network configurations than expected by chance. In turn, this is compatible with a limited number of states as 
revealed by a HMM analysis (Fig. 7). 

Note that a HMM analysis on 30 simulated neurons taken from 30 different clusters is bound to produce a 
larger number of states than detected in the data, because in the latter only ensembles of 3 to 9 simultaneously 
recorded neurons were available. Moreover, it is not guaranteed that the recorded neurons all came from 
different clusters. Thus, to facilitate comparison between model and data, we picked a distribution of ensemble 
sizes matching the empirical distribution (3 to 9 neurons), and chose standard physiological values for all other 
parameters (Appendix 2). Three representative trials from an ensemble with 9 neurons are shown in Fig. 8A, 
revealing the re-occurrence of HMM states across trials, with 8.3 ±1.7 states (mean±SD), ranging from 5 to 
12 (a range much closer to the range of 3 to 7 states found in the data), and approximately exponential duration 
with mean of 239 ms (95% Cl: [234, 243]). The number of states detected was correlated with the number of 
neurons in each state (i? 2 = 0.36, p < 0.01), as could be expected since the probability of detecting a state 
will increase in larger ensembles until neurons from all clusters are sampled. The model states matched other 
characteristic features of the data: state transitions followed CPs with a similar shape of the CPTA (Fig. 8B; 
compare with Fig. 3B) and were often congruent with the co-modulation of a fraction of the neurons (Fig. 8C; 
compare with Fig. 3C). Moreover, we found that 44% of the neurons had 3 or more firing rates (Fig. 8D), 
essentially the same fraction (42%) found in the data (see Fig. 3F). No additional tuning of the parameters was 
required to obtain this match. 

3.6 Comparison between ongoing and stimulus-evoked activity: reduction of multi-stability 

Previous work has proposed various relationships between patterns of ongoing activity and responses evoked 
by sensory stimuli, see e.g. [18, 20, 24, 67-70]. We performed a series of analyses of the experimental and 
simulated data to determine whether the model network, developed to reproduce GC ongoing activity, captured 
the essential features of taste-evoked activity with no additional tuning of parameters. We first performed an 
HMM analysis on the data recorded during evoked activity in the [0, 5] s interval post-taste delivery, as in 
[8], see Fig. 9A. We found a range of 4 to 11 states per taste across sessions (mean±SD: 7.2 ± 1.6) with an 
approximate exponential distribution of durations with mean 306 ms (95% Cl: [278, 334]). However, during 
evoked activity only 8% of the neurons had more than 2 different firing rates across states compared to 42% 
during ongoing activity (Fig. 9B, % 2 (1) = 51, p < 10~ 12 ). This suggests that a stimulus steers the firing 
rates of single neurons away from the multi-stable regime, hence reducing the range and number of different 
firing rates available. Moreover, in keeping with previous reports [42], the stimulus caused a significant drop in 
trial-by-trial variability, as measured by the mean-matched Fano factor (see Appendix 2 and Fig. 9C). 

To gain insight into the mechanism responsible for the observed reduction in multi-stability, we investigated 
whether this new configuration imposed by the stimulus is compatible with the attractor landscape predicted by 
the model (Fig. 4B). We analyzed the behavior of the model network in the presence of a stimulus mimicking 
thalamic input following taste delivery (Appendix 2). Four stimuli were used (to mimic the 4 tastants used 
in experiment), which differed only in the neurons they targeted, according to the same empirical distribution 
found in the data (2 representative stimuli are shown in the cartoon of Fig. 10A). Specifically, the fractions of 
neurons responsive to n = 1, 2, 3 or all 4 stimuli were 17% (27/162), 22% (36/162), 26% (42/162), and 35% 
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Figure 8. Ongoing activity in the model is characterized by sequences of states. A: three representative trials of simulated 
data after keeping only 9 randomly chosen neurons, together with a representation of the corresponding HMM states (lower 
panel); notations as in panel A. B: CP-triggered average (CPTA) of ensemble transitions for the simulated dataset (black) 
and for a trial-shuffled dataset (blue); same conventions as in Fig. 3B. Thick black line: p < 0.05 (t-test with multiple-bin 
Bonferroni correction). X-axis: time lag from CP (s); Y-axis: probability of a transition given a CP at t = 0. C: fraction 
of transitions co-occuring with CPs within the 200 ms significant window from panel B. X-axis: number of single neurons 
CPs divided by the ensemble size; Y-axis: fraction of transitions co-occurring with CPs (%). D: number of different firing 
rates per neuron across all sessions and neurons: 44% of all model neurons had three or more different firing rates across 
ensemble states (compare with Fig. 3F). X-axis: minimal number of different firing rates across states; Y-axis: fraction of 
neurons (%). 


(57/162), respectively. All other model parameters were chosen as in the analysis of ongoing activity. With 
a biologically realistic stimulus peaking at 30% of v exl (see Appendix 2), we found a range of 4 to 17 states 
per stimulus across sessions (meaniSD: 10.0 ± 3.0), with an approximate exponential distribution of state 
durations with mean 227 ms (95% Cl: [219,235] ms). Representative trials from simulated activity evoked by 
two of the four different stimuli are shown in Fig. 10A, together with their segmentation in HMM states. As 
in the data, a significantly smaller fraction (15%) of neurons had 3 or more different firing rates across states 
compared to ongoing activity (44%; Fig. 10B, x 2 (l) = 24, p < 10 -5 ). This result was robust to changes in 
stimulus shape and amplitude; using different stimuli gave a comparable fraction of multi-stable neurons (either 
varying stimulus peak amplitude from 10% to 30% of zz ext , or using a box stimulus; not shown, see Appendix 2). 
The model also captured the stimulus-induced reduction in trial-by-trial variability found in the data (Fig. 10C), 
as measured by the mean-matched Fano Factor, confirming previous results found with similar model networks 
[13, 15]. 

Both the empirical data and the behavior of the network under stimulation could be explained with the 
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Figure 9. Taste-evoked activity: reduction of multi-stability and trial-by-trial variability in the experimental data. A: two 
representative examples of population rasters together with their segmentation in HMM states during evoked activity in 
GC, in response to citric acid (left) and sucrose (right) delivery, respectively (the red arrow marks taste delivery at t = 0; 
remaining notations as in Fig. 1A). Upper panel: X-axis, time relative to taste delivery (0 = stimulus delivery); X-axis, 
left, ensemble neuron index, right, posterior probability of FIMM states. Lower panel: X-axis for lower panel: firing rates 
(spks/s); X-axis, ensemble neuron index. B: fraction of multi-stable neurons across all states during ongoing (left) and 
evoked activity (right) in GC (** = p < 0.01, \ 2 test). A'-axis: ongoing and evoked conditions; Y-axis: fraction of 
multi-stable neurons (%). C: time course of the mean-matched Fano Factor (FF) in a time interval around taste delivery 
(occurring at time f = 0) across all data. Shaded bounds: 95% CIs. The thick horizontal black line indicates bins where the 
evoked FF is significantly different from baseline (p < 0.05, see Appendix 2). X-axis: time (s); X-axis: Fano Factor. 


same theoretical analysis used to predict ongoing activity (Fig. 4). We repeated the mean field analysis for 
a fixed value of J + = 5.2 (green vertical line and arrow in Fig. 4B) but varying the stimulus amplitude. 
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Figure 10. Taste-evoked activity: reduction of multi-stability and trial-by trial variability in the model. A: two examples of 
population rasters together with their segmentation in HMM states during evoked activity in the model network, in response 
to two surrogate stimuli. Clusters may be selective to more than one stimulus (central lower panel: dark circles and colored 
areas represent, respectively, individual clusters and their selectivity to the two different stimuli A (blue) and B (red)); 
remaining conventions as in Fig. 9A. Model stimuli (central upper panel) consist of double exponentials peaking at 30% of 
Vext (50 ms rise time and 500 ms decay time). X-axis: time (s), X-axis, stimulus amplitude (%). B: fraction of multi-stable 
neurons across all states during ongoing and evoked activity in the model (** = p < 0.01). A'-axis: ongoing and evoked 
conditions; X-axis: fraction of multi-stable neurons (%). C: time course of the mean-matched Fano Factor in a time interval 
around stimulus presentation (occurring at time t = 0) across all simulated sessions (same conventions as in Fig. 9C). 
X-axis: time (s); X-axis: Fano factor. 


The results are shown in Fig. 11 A. The number of active clusters tends to increase with stimulus amplitude. 
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Figure 11. Mean field predictions of evoked activity and comparison with model simulations. A: attractor states of the model 
during evoked activity predicted by mean field theory for J+ = 5.2 and as a function of stimulus amplitude. Diamonds: 
mean firing rate in active clusters for each attractor state, featuring the number of active clusters indicated by the red 
numbers (see also Fig. 4B). Stimulus amplitude varied from 1% to 45% of u ext . Vertical dashed lines: attractor states for 
stimuli of 5%, 25% and 45%. The firing rate range decreases as the stimulus amplitude is increased, and it reduces to a 
single configuration with 15 active clusters for stimuli > 30%. X-axis: stimulus increase (%); X-axis: firing rate (spks/s). 
B: mean firing rates in active clusters in simulations of the model for the three stimuli marked by vertical lines in panel 
A (vertical bars: S.D.). A'-axis: number of active clusters; X-axis: average cluster firing rate (spks/s). Inset: frequency 
of configurations with different numbers of active clusters for 5% stimulus. X-axis: active clusters; X-axis: frequency 
of occurrence (% of total simulation time across all simulations). C: rasterplot of spike trains from excitatory clustered 
neurons (background excitatory population not shown) in a representative trial encompassing both periods of ongoing and 
evoked activity (compare to Fig. 6A), using a box stimulus with amplitude at to 30% of v ex t in the [0, 2] s interval (ongoing 
activity corresponds to the [—2,0] s interval). Vertical red line: stimulus onset. X-axis: time (s); X-axis: neuron index. 
D: number of active clusters vs. time (representative trial, panel C; same notation as in Fig. 6B). X-axis: time (s); X-axis: 
active clusters. 
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while the configurations with fewer active clusters gradually disappear. Moreover, one observes a reduction in 
the range of firing rates across states until, for a stimulus amplitude > 30% of only configurations with 
15 active clusters are theoretically possible. Computer simulations confirmed these predictions (Fig. 1 IB): 
although the match between predicted and simulated firing rates was not perfect (presumably due to finite size 
effects and the approximations used for synaptic dynamics), the predicted narrowing of the firing rate range 
for stronger stimuli was evident. As with ongoing activity, we also found that the network spent most of 
its time visiting configurations with an intermediate number of active clusters among all those theoretically 
possible predicted in mean field (inset of Fig. 1 IB), suggesting, again, that different network configurations 
have different basins of attraction, with larger basins corresponding to longer durations. An example simulation 
for a stimulus amplitude of 30% of // ex , is shown in Fig. 1 1C, together with the running score of the number of 
active clusters (Fig. 11D). In conclusion, stimulating a cortical network during a state of multi-stable ongoing 
activity reduces not only the trial-by-trial variability, as previously reported, but also multi-stability, defined 
as the single neurons ability to exhibit multiple (3 or more) firing rates across states. This may be linked to a 
reduction in complexity of evoked neural representations [71], 

4 Discussion 

Taste-evoked activity in GC has been characterized as a progression through a sequence of metastable states, 
where each state is a collection of firing rates across simultaneously recorded neurons [8]. Here we demonstrate 
that metastable states can also be observed during ongoing activity and report several quantitative comparisons 
between ongoing and evoked activity. The most notable difference is that single neurons are multi-stable during 
ongoing activity (expressing 3 or more different firing rates across states), whereas they are mostly bi-stable 
during evoked activity (expressing at most 2 different firing rates across states). These results were reproduced 
in a biologically plausible, spiking network model of cortical activity amenable to theoretical analysis. The 
network had dense and sufficiently strong intra-cluster connections, which endowed it with a rich landscape of 
stable states with several different firing rates. To our knowledge, this rich repertoire of network attractors and 
its specific modification under stimulation has not been reported before. The stable states become metastable in 
a finite network, wherein endogenously generated fluctuations induce an itinerant dynamics among ensemble 
states. As observed in the data, this phenomenon did not require overt external stimulation. The model captured 
other essential features of the data, such as the distribution of single neurons firing rates across states and the 
partial co-activation of neurons associated to each ensemble transition. Moreover, without additional tuning, the 
model reproduced the reduction of multi-stability and trial-by-trial variability at stimulus onset. The reduction 
in trial-by-trial variability confirms previous results obtained across several brain areas [42], which were also 
reproduced in spiking models similar to ours [13, 15, 17]. As for multi-stability, our theory offers a mechanistic 
explanation of its origin and its stimulus-induced reduction. The theory predicts the complete abolition of 
multi-stability above a critical value of the external stimulation (Fig. 11 A), when the network exhibits only one 
firing rate across all neurons in active clusters. In the more plausible range of mid-strength stimulations, the 
model shows a reduction in multi-stability and a contraction of the range of firing rates. 

4.1 External and internal sources of metastability in GC 

Metastable activity can be induced by an external stimulus as recently shown by [10, 72] in a spiking network 
model. Their model relies on feed-forward connections between recurrent cortical modules responsible for 
taste detection, taste identification, and decision to spit or swallow. Each stage corresponds to the activation 
of a particular module, and transitions between stages are driven by stochastic fluctuations in the network. 
While this model reproduces the dynamics of taste-evoked activity, it relies on an external stimulus to ignite 
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the transitions, and thus would not explain metastable dynamics during ongoing activity. Our model provides 
a mechanistic explanation of intrinsically generated state sequences via a common mechanism during both 
ongoing and evoked activity. The segmentation of ongoing activity in a sequence of metastable states suggests 
that the network hops among persistent activity configurations that lose stability by mechanisms including 
intrinsically generated noise [13-17, 73], some process of adaptation [64, 73-75], or winnerless competition 
[76]. In our model and in previous models similar to ours [13, 15, 17], the mechanism is internally generated 
noise that destabilizes the stable attractors predicted by mean field theory. In a finite network, the stable states 
become metastable due to fluctuations of neural activity similar to those typically observed in vivo [56, 57, 77]. 
These fluctuations are a consequence of the network being in the balanced regime [50, 51, 65], which produces 
irregular and asynchronous spiking activity [45, 52], In the network containing a finite number of neurons, 
this fluctuating activity can destabilize the stable states and produce hopping behavior among metastable states 
[13-17], 

4.2 Relationship between ongoing and evoked activity 

The nature of the relationship between ongoing and evoked states remains elusive. One proposal for the role 
of ongoing activity is that it serves as a repertoire of representations sampled by the neural circuit during 
evoked activity [18, 20, 24], According to this proposal, evoked states occupy a subset of ongoing states 
in a reduced representational space, such as the space of principal components, or the space obtained after 
multidimensional scaling [24]. This applies especially to the activity of auditory or somatosensory neurons 
during high and low activity states called “UP” and “DOWN” states [24, 78-80]. An alternative proposal 
regards ongoing activity as a Bayesian prior that incrementally adjusts to the statistics of external stimuli during 
early development [21, 70]. In this framework, ongoing activity starts out as initially different from evoked 
activity, and progressively shifts towards it. Our models explanation of the genesis and structure of ensemble 
states has two main implications: if, on the one hand, ongoing and evoked activities are undoubtedly linked 
(they emerge from the same network structure and organization); on the other hand they may differ in important 
ways. Given the networks structure, only a handful of states are allowed to emerge through network dynamics, 
which might suggest that ongoing activity constrains the repertoire of sensory responses. However, the same 
networks structure, especially its synaptic organization in potentiated clusters, can be obtained as the result 
of learning external stimuli [17, 43, 81, 82], so from this point of view it is the realm of evoked responses 
that defines the structure of ongoing activity. In any case, even though originating from the same network 
structure and mechanisms, evoked and ongoing activity can be characteristically different. They may differ 
in the number of ensemble states or firing rate distributions across neurons, resulting in only a partial overlap 
when analyzed as in [24] (not shown). Stimuli drive neurons at higher and lower firing rates than typically 
observed in ongoing activity, preventing ongoing activity from completely encompassing the evoked states. 
This is a direct consequence of our finding that activity in GC is not compatible with UP and DOWN states, 
both in terms of the dynamics of state transitions (Fig. 3D) and in terms of the wider range of firing rates 
attained by single neurons (multi-stability. Fig. 3F). The observed differences between ongoing and evoked 
activity can be understood thanks to our spiking network model. This model shows how the presence of a 
stimulus induces a change in the landscape of metastable configurations (Fig. 11 A). This change depends on 
the strength of the stimulus and is incremental, quenching the range and the number of the active clusters in 
admissible configurations. Beyond a critical point (> 30% stimulus strength in Fig. 1 1A), the only states left 
are in configurations where all stimulated clusters are simultaneously active, and evoked states are intrinsically 
different from ongoing states. Although they still relate to ongoing states due to their common network origin, 
evoked states contain information that is not available in the ongoing activity itself (and indeed, this information 
can be used to decode the taste, see e.g. [8, 27, 83]). 
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4.3 Mechanisms of itinerant multi-stable activity 

We have shown that our network model, although balanced, produces coordinated co-activation across clusters, 
mostly because of recurrent inhibition and heterogeneity in cluster size. Recurrent inhibition reduces the num¬ 
ber of clusters that can be simultaneously active, whereas heterogeneity in cluster size inflates the probability of 
visiting larger clusters at the expense of smaller ones, thus initiating recurrent sequences of states. In turn, this 
results in more coordinated cluster dynamics and a much more limited number of network configurations than 
expected by chance. It is tempting to infer, on the basis of the results presented in Figures 3 and 8, a similar 
degree of partial coordination in the experimental data during ongoing activity. Given the limited size of our en¬ 
sembles and the limited number of neurons participating in state transitions, it is difficult to determine whether 
the experimental data reflect the same degree of coordination predicted by the model. While recent evidence 
of clustered subpopulations in the frontal and motor cortex of behaving monkeys seem to support our model 
[84], future experiments making use of high-density recordings may provide more accurate measurements of 
the degree of coordination among clusters. It is worth noting that the latter may not be fixed but rather depend 
on the behavioral state of the animal. During sleep or anesthesia, a spiking network may experience more syn¬ 
chronous global states, while during alertness it may exhibit intermediate degrees of coordination such as those 
predicted by our model. Similarly, evoked activity might display a more global degree of coordination than 
ongoing activity. Depending on the strength of inter-cluster connections and the size of individual clusters, our 
model can produce networks with different degrees of coordination and thus provide a link between behavioral 
states and mechanisms of coordination. Finally, our model offers a novel mechanism for multi-stable firing 
rates across states, i.e., the observation that approximately half of our neurons have 3 or more firing rates across 
states. Even though a network with heterogeneous synaptic weights could have a distribution of different firing 
rates across different neurons [85], each neuron would only be firing at two different firing rates (bi-stability). 
This is at odds with multi-stability. In our model, having several firing rate patterns in each single neuron im¬ 
plies that the network visits states with different numbers of active clusters, where the firing rate of each cluster 
depends on the number of active clusters at any given time. 

State sequences seem to encode gustatory information and are believed to play a role in taste processing 
and taste-guided decisions [8, 9, 86], Potential mechanisms to read out such sequences in populations of 
spiking neurons are being investigated [87]. In addition to taste coding, multi-stable activity across a variety of 
metastable states may enrich the network’s ability to encode high dimensional and temporally dynamic sensory 
experiences. The investigation of this possibility and its link to the dimensionality of neural representations 
[71] is left for future work. 
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